home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 21 / AACD 21.iso / AACD / Utilities / Ghostscript / src / gsfemu.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-01-01  |  19.7 KB  |  851 lines
  1. /* Copyright (C) 1989, 1996, 1998 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of AFPL Ghostscript.
  4.   
  5.   AFPL Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author or
  6.   distributor accepts any responsibility for the consequences of using it, or
  7.   for whether it serves any particular purpose or works at all, unless he or
  8.   she says so in writing.  Refer to the Aladdin Free Public License (the
  9.   "License") for full details.
  10.   
  11.   Every copy of AFPL Ghostscript must include a copy of the License, normally
  12.   in a plain ASCII text file named PUBLIC.  The License grants you the right
  13.   to copy, modify and redistribute AFPL Ghostscript, but only under certain
  14.   conditions described in the License.  Among other things, the License
  15.   requires that the copyright notice and this notice be preserved on all
  16.   copies.
  17. */
  18.  
  19. /*$Id: gsfemu.c,v 1.2 2000/09/19 19:00:28 lpd Exp $ */
  20. /* Floating point emulator for gcc */
  21.  
  22. /* We actually only need arch.h + uint and ulong, but because signal.h */
  23. /* may include <sys/types.h>, we have to include std.h to handle (avoid) */
  24. /* redefinition of type names. */
  25. #include "std.h"
  26. #include <signal.h>
  27.  
  28. /*#define TEST */
  29.  
  30. /*
  31.  * This package is not a fully general IEEE floating point implementation.
  32.  * It implements only one rounding mode (round to nearest) and does not
  33.  * generate or properly handle NaNs.
  34.  *
  35.  * The names and interfaces of the routines in this package were designed to
  36.  * work with gcc.  This explains some peculiarities such as the presence of
  37.  * single-precision arithmetic (forbidden by the ANSI standard) and the
  38.  * omission of unsigned-to-float conversions (which gcc implements with
  39.  * truly grotesque inline code).
  40.  *
  41.  * The following routines do not yet handle denormalized numbers:
  42.  *      addd3 (operands or result)
  43.  *      __muldf3 (operands or result)
  44.  *      __divdf3 (operands or result)
  45.  *      __truncdfsf2 (result)
  46.  *      __extendsfdf2 (operand)
  47.  */
  48.  
  49. /*
  50.  * IEEE single-precision floats have the format:
  51.  *      sign(1) exponent(8) mantissa(23)
  52.  * The exponent is biased by +127.
  53.  * The mantissa is a fraction with an implicit integer part of 1,
  54.  * unless the exponent is zero.
  55.  */
  56. #define fx(ia) (((ia) >> 23) & 0xff)
  57. #define fx_bias 127
  58. /*
  59.  * IEEE double-precision floats have the format:
  60.  *      sign(1) exponent(11) mantissa(52)
  61.  * The exponent is biased by +1023.
  62.  */
  63. #define dx(ld) ((ld[msw] >> 20) & 0x7ff)
  64. #define dx_bias 1023
  65.  
  66. #if arch_is_big_endian
  67. #  define msw 0
  68. #  define lsw 1
  69. #else
  70. #  define msw 1
  71. #  define lsw 0
  72. #endif
  73. /* Double arguments/results */
  74. #define la ((const long *)&a)
  75. #define ula ((const ulong *)&a)
  76. #define lb ((const long *)&b)
  77. #define ulb ((const ulong *)&b)
  78. #define dc (*(const double *)lc)
  79. /* Float arguments/results */
  80. #define ia (*(const long *)&a)
  81. #define ua (*(const ulong *)&a)
  82. #define ib (*(const long *)&b)
  83. #define ub (*(const ulong *)&b)
  84. #define fc (*(const float *)&lc)
  85.  
  86. /* Round a double-length fraction by adding 1 in the lowest bit and */
  87. /* then shifting right by 1. */
  88. #define roundr1(ms, uls)\
  89.   if ( uls == 0xffffffff ) ms++, uls = 0;\
  90.   else uls++;\
  91.   uls = (uls >> 1) + (ms << 31);\
  92.   ms >>= 1
  93.  
  94. /* Extend a non-zero float to a double. */
  95. #define extend(lc, ia)\
  96.   ((lc)[msw] = ((ia) & 0x80000000) + (((ia) & 0x7fffffff) >> 3) + 0x38000000,\
  97.    (lc)[lsw] = (ia) << 29)
  98.  
  99. /* ---------------- Arithmetic ---------------- */
  100.  
  101. /* -------- Add/subtract/negate -------- */
  102.  
  103. double
  104. __negdf2(double a)
  105. {
  106.     long lc[2];
  107.  
  108.     lc[msw] = la[msw] ^ 0x80000000;
  109.     lc[lsw] = la[lsw];
  110.     return dc;
  111. }
  112. float
  113. __negsf2(float a)
  114. {
  115.     long lc = ia ^ 0x80000000;
  116.  
  117.     return fc;
  118. }
  119.  
  120. double
  121. __adddf3(double a, double b)
  122. {
  123.     long lc[2];
  124.     int expt = dx(la);
  125.     int shift = expt - dx(lb);
  126.     long sign;
  127.     ulong msa, lsa;
  128.     ulong msb, lsb;
  129.  
  130.     if (shift < 0) {        /* Swap a and b so that expt(a) >= expt(b). */
  131.     double temp = a;
  132.  
  133.     a = b, b = temp;
  134.     expt += (shift = -shift);
  135.     }
  136.     if (shift >= 54)        /* also picks up most cases where b == 0 */
  137.     return a;
  138.     if (!(lb[msw] & 0x7fffffff))
  139.     return a;
  140.     sign = la[msw] & 0x80000000;
  141.     msa = (la[msw] & 0xfffff) + 0x100000, lsa = la[lsw];
  142.     msb = (lb[msw] & 0xfffff) + 0x100000, lsb = lb[lsw];
  143.     if ((la[msw] ^ lb[msw]) >= 0) {    /* Adding numbers of the same sign. */
  144.     if (shift >= 32)
  145.         lsb = msb, msb = 0, shift -= 32;
  146.     if (shift) {
  147.         --shift;
  148.         lsb = (lsb >> shift) + (msb << (32 - shift));
  149.         msb >>= shift;
  150.         roundr1(msb, lsb);
  151.     }
  152.     if (lsb > (ulong) 0xffffffff - lsa)
  153.         msa++;
  154.     lsa += lsb;
  155.     msa += msb;
  156.     if (msa > 0x1fffff) {
  157.         roundr1(msa, lsa);
  158.         /* In principle, we have to worry about exponent */
  159.         /* overflow here, but we don't. */
  160.         ++expt;
  161.     }
  162.     } else {            /* Adding numbers of different signs. */
  163.     if (shift > 53)
  164.         return a;        /* b can't affect the result, even rounded */
  165.     if (shift == 0 && (msb > msa || (msb == msa && lsb >= lsa))) {    /* This is the only case where the sign of the result */
  166.         /* differs from the sign of the first operand. */
  167.         sign ^= 0x80000000;
  168.         msa = msb - msa;
  169.         if (lsb < lsa)
  170.         msa--;
  171.         lsa = lsb - lsa;
  172.     } else {
  173.         if (shift >= 33) {
  174.         lsb = ((msb >> (shift - 32)) + 1) >> 1;        /* round */
  175.         msb = 0;
  176.         } else if (shift) {
  177.         lsb = (lsb >> (shift - 1)) + (msb << (33 - shift));
  178.         msb >>= shift - 1;
  179.         roundr1(msb, lsb);
  180.         }
  181.         msa -= msb;
  182.         if (lsb > lsa)
  183.         msa--;
  184.         lsa -= lsb;
  185.     }
  186.     /* Now renormalize the result. */
  187.     /* For the moment, we do this the slow way. */
  188.     if (!(msa | lsa))
  189.         return 0;
  190.     while (msa < 0x100000) {
  191.         msa = (msa << 1) + (lsa >> 31);
  192.         lsa <<= 1;
  193.         expt -= 1;
  194.     }
  195.     if (expt <= 0) {    /* Underflow.  Return 0 rather than a denorm. */
  196.         lc[msw] = sign;
  197.         lc[lsw] = 0;
  198.         return dc;
  199.     }
  200.     }
  201.     lc[msw] = sign + ((ulong) expt << 20) + (msa & 0xfffff);
  202.     lc[lsw] = lsa;
  203.     return dc;
  204. }
  205. double
  206. __subdf3(double a, double b)
  207. {
  208.     long nb[2];
  209.  
  210.     nb[msw] = lb[msw] ^ 0x80000000;
  211.     nb[lsw] = lb[lsw];
  212.     return a + *(const double *)nb;
  213. }
  214.  
  215. float
  216. __addsf3(float a, float b)
  217. {
  218.     long lc;
  219.     int expt = fx(ia);
  220.     int shift = expt - fx(ib);
  221.     long sign;
  222.     ulong ma, mb;
  223.  
  224.     if (shift < 0) {        /* Swap a and b so that expt(a) >= expt(b). */
  225.     long temp = ia;
  226.  
  227.     *(long *)&a = ib;
  228.     *(long *)&b = temp;
  229.     expt += (shift = -shift);
  230.     }
  231.     if (shift >= 25)        /* also picks up most cases where b == 0 */
  232.     return a;
  233.     if (!(ib & 0x7fffffff))
  234.     return a;
  235.     sign = ia & 0x80000000;
  236.     ma = (ia & 0x7fffff) + 0x800000;
  237.     mb = (ib & 0x7fffff) + 0x800000;
  238.     if ((ia ^ ib) >= 0) {    /* Adding numbers of the same sign. */
  239.     if (shift) {
  240.         --shift;
  241.         mb = ((mb >> shift) + 1) >> 1;
  242.     }
  243.     ma += mb;
  244.     if (ma > 0xffffff) {
  245.         ma = (ma + 1) >> 1;
  246.         /* In principle, we have to worry about exponent */
  247.         /* overflow here, but we don't. */
  248.         ++expt;
  249.     }
  250.     } else {            /* Adding numbers of different signs. */
  251.     if (shift > 24)
  252.         return a;        /* b can't affect the result, even rounded */
  253.     if (shift == 0 && mb >= ma) {
  254.         /* This is the only case where the sign of the result */
  255.         /* differs from the sign of the first operand. */
  256.         sign ^= 0x80000000;
  257.         ma = mb - ma;
  258.     } else {
  259.         if (shift) {
  260.         --shift;
  261.         mb = ((mb >> shift) + 1) >> 1;
  262.         }
  263.         ma -= mb;
  264.     }
  265.     /* Now renormalize the result. */
  266.     /* For the moment, we do this the slow way. */
  267.     if (!ma)
  268.         return 0;
  269.     while (ma < 0x800000) {
  270.         ma <<= 1;
  271.         expt -= 1;
  272.     }
  273.     if (expt <= 0) {
  274.         /* Underflow.  Return 0 rather than a denorm. */
  275.         lc = sign;
  276.         return fc;
  277.     }
  278.     }
  279.     lc = sign + ((ulong)expt << 23) + (ma & 0x7fffff);
  280.     return fc;
  281. }
  282.  
  283. float
  284. __subsf3(float a, float b)
  285. {
  286.     long lc = ib ^ 0x80000000;
  287.  
  288.     return a + fc;
  289. }
  290.  
  291. /* -------- Multiplication -------- */
  292.  
  293. double
  294. __muldf3(double a, double b)
  295. {
  296.     long lc[2];
  297.     ulong sign;
  298.     uint H, I, h, i;
  299.     ulong p0, p1, p2;
  300.     ulong expt;
  301.  
  302.     if (!(la[msw] & 0x7fffffff) || !(lb[msw] & 0x7fffffff))
  303.     return 0;
  304.     /*
  305.      * We now have to multiply two 53-bit fractions to produce a 53-bit
  306.      * result.  Since the idiots who specified the standard C libraries
  307.      * don't allow us to use the 32 x 32 => 64 multiply that every
  308.      * 32-bit CPU provides, we have to chop these 53-bit numbers up into
  309.      * 14-bit chunks so we can use 32 x 32 => 32 multiplies.
  310.      */
  311. #define chop_ms(ulx, h, i)\
  312.   h = ((ulx[msw] >> 7) & 0x1fff) | 0x2000,\
  313.   i = ((ulx[msw] & 0x7f) << 7) | (ulx[lsw] >> 25)
  314. #define chop_ls(ulx, j, k)\
  315.   j = (ulx[lsw] >> 11) & 0x3fff,\
  316.   k = (ulx[lsw] & 0x7ff) << 3
  317.     chop_ms(ula, H, I);
  318.     chop_ms(ulb, h, i);
  319. #undef chop
  320. #define prod(m, n) ((ulong)(m) * (n))    /* force long multiply */
  321.     p0 = prod(H, h);
  322.     p1 = prod(H, i) + prod(I, h);
  323.     /* If these doubles were produced from floats or ints, */
  324.     /* all the other products will be zero.  It's worth a check. */
  325.     if ((ula[lsw] | ulb[lsw]) & 0x1ffffff) {    /* No luck. */
  326.     /* We can almost always avoid computing p5 and p6, */
  327.     /* but right now it's not worth the trouble to check. */
  328.     uint J, K, j, k;
  329.  
  330.     chop_ls(ula, J, K);
  331.     chop_ls(ulb, j, k);
  332.     {
  333.         ulong p6 = prod(K, k);
  334.         ulong p5 = prod(J, k) + prod(K, j) + (p6 >> 14);
  335.         ulong p4 = prod(I, k) + prod(J, j) + prod(K, i) + (p5 >> 14);
  336.         ulong p3 = prod(H, k) + prod(I, j) + prod(J, i) + prod(K, h) +
  337.         (p4 >> 14);
  338.  
  339.         p2 = prod(H, j) + prod(I, i) + prod(J, h) + (p3 >> 14);
  340.     }
  341.     } else {
  342.     p2 = prod(I, i);
  343.     }
  344.     /* Now p0 through p2 hold the result. */
  345. /****** ASSUME expt IS 32 BITS WIDE. ******/
  346.     expt = (la[msw] & 0x7ff00000) + (lb[msw] & 0x7ff00000) -
  347.     (dx_bias << 20);
  348.     /* Now expt is in the range [-1023..3071] << 20. */
  349.     /* Values outside the range [1..2046] are invalid. */
  350.     p1 += p2 >> 14;
  351.     p0 += p1 >> 14;
  352.     /* Now the 56-bit result consists of p0, the low 14 bits of p1, */
  353.     /* and the low 14 bits of p2. */
  354.     /* p0 can be anywhere between 2^26 and almost 2^28, so we might */
  355.     /* need to adjust the exponent by 1. */
  356.     if (p0 < 0x8000000) {
  357.     p0 = (p0 << 1) + ((p1 >> 13) & 1);
  358.     p1 = (p1 << 1) + ((p2 >> 13) & 1);
  359.     p2 <<= 1;
  360.     } else
  361.     expt += 0x100000;
  362.     /* The range of expt is now [-1023..3072] << 20. */
  363.     /* Round the mantissa to 53 bits. */
  364.     if (!((p2 += 4) & 0x3ffc) && !(++p1 & 0x3fff) && ++p0 >= 0x10000000) {
  365.     p0 >>= 1, p1 = 0x2000;
  366.     /* Check for exponent overflow, just in case. */
  367.     if ((ulong) expt < 0xc0000000)
  368.         expt += 0x100000;
  369.     }
  370.     sign = (la[msw] ^ lb[msw]) & 0x80000000;
  371.     if (expt == 0) {        /* Underflow.  Return 0 rather than a denorm. */
  372.     lc[msw] = sign;
  373.     lc[lsw] = 0;
  374.     } else if ((ulong) expt >= 0x7ff00000) {    /* Overflow or underflow. */
  375.     if ((ulong) expt <= 0xc0000000) {    /* Overflow. */
  376.         raise(SIGFPE);
  377.         lc[msw] = sign + 0x7ff00000;
  378.         lc[lsw] = 0;
  379.     } else {        /* Underflow.  See above. */
  380.         lc[msw] = sign;
  381.         lc[lsw] = 0;
  382.     }
  383.     } else {
  384.     lc[msw] = sign + expt + ((p0 >> 7) & 0xfffff);
  385.     lc[lsw] = (p0 << 25) | ((p1 & 0x3fff) << 11) | ((p2 >> 3) & 0x7ff);
  386.     }
  387.     return dc;
  388. #undef prod
  389. }
  390. float
  391. __mulsf3(float a, float b)
  392. {
  393.     uint au, al, bu, bl, cu, cl, sign;
  394.     long lc;
  395.     uint expt;
  396.  
  397.     if (!(ia & 0x7fffffff) || !(ib & 0x7fffffff))
  398.     return 0;
  399.     au = ((ia >> 8) & 0x7fff) | 0x8000;
  400.     bu = ((ib >> 8) & 0x7fff) | 0x8000;
  401.     /* Since 0x8000 <= au,bu <= 0xffff, 0x40000000 <= cu <= 0xfffe0001. */
  402.     cu = au * bu;
  403.     if ((al = ia & 0xff) != 0) {
  404.     cl = bu * al;
  405.     } else
  406.     cl = 0;
  407.     if ((bl = ib & 0xff) != 0) {
  408.     cl += au * bl;
  409.     if (al)
  410.         cl += (al * bl) >> 8;
  411.     }
  412.     cu += cl >> 8;
  413.     sign = (ia ^ ib) & 0x80000000;
  414.     expt = (ia & 0x7f800000) + (ib & 0x7f800000) - (fx_bias << 23);
  415.     /* expt is now in the range [-127..383] << 23. */
  416.     /* Values outside [1..254] are invalid. */
  417.     if (cu <= 0x7fffffff)
  418.     cu <<= 1;
  419.     else
  420.     expt += 1 << 23;
  421.     cu = ((cu >> 7) + 1) >> 1;
  422.     if (expt < 1 << 23)
  423.     lc = sign;        /* underflow */
  424.     else if (expt > (uint)(254 << 23)) {
  425.     if (expt <= 0xc0000000) { /* overflow */
  426.         raise(SIGFPE);
  427.         lc = sign + 0x7f800000;
  428.     } else {        /* underflow */
  429.         lc = sign;
  430.     }
  431.     } else
  432.     lc = sign + expt + cu - 0x800000;
  433.     return fc;
  434. }
  435.  
  436. /* -------- Division -------- */
  437.  
  438. double
  439. __divdf3(double a, double b)
  440. {
  441.     long lc[2];
  442.  
  443.     /*
  444.      * Multi-precision division is really, really awful.
  445.      * We generate the result more or less brute force,
  446.      * 11 bits at a time.
  447.      */
  448.     ulong sign = (la[msw] ^ lb[msw]) & 0x80000000;
  449.     ulong msa = (la[msw] & 0xfffff) | 0x100000, lsa = la[lsw];
  450.     ulong msb = (lb[msw] & 0xfffff) | 0x100000, lsb = lb[lsw];
  451.     uint qn[5];
  452.     int i;
  453.     ulong msq, lsq;
  454.     int expt = dx(la) - dx(lb) + dx_bias;
  455.  
  456.     if (!(lb[msw] & 0x7fffffff)) {    /* Division by zero. */
  457.     raise(SIGFPE);
  458.     lc[lsw] = 0;
  459.     lc[msw] =
  460.         (la[msw] & 0x7fffffff ?
  461.          sign + 0x7ff00000 /* infinity */ :
  462.          0x7ff80000 /* NaN */ );
  463.     return dc;
  464.     }
  465.     if (!(la[msw] & 0x7fffffff))
  466.     return 0;
  467.     for (i = 0; i < 5; ++i) {
  468.     uint q;
  469.     ulong msp, lsp;
  470.  
  471.     msa = (msa << 11) + (lsa >> 21);
  472.     lsa <<= 11;
  473.     q = msa / msb;
  474.     /* We know that 2^20 <= msb < 2^21; q could be too high by 1, */
  475.     /* but it can't be too low. */
  476.     /* Set p = d * q. */
  477.     msp = q * msb;
  478.     lsp = q * (lsb & 0x1fffff);
  479.     {
  480.         ulong midp = q * (lsb >> 21);
  481.  
  482.         msp += (midp + (lsp >> 21)) >> 11;
  483.         lsp += midp << 21;
  484.     }
  485.     /* Set a = a - p, i.e., the tentative remainder. */
  486.     if (msp > msa || (lsp > lsa && msp == msa)) {    /* q was too large. */
  487.         --q;
  488.         if (lsb > lsp)
  489.         msp--;
  490.         lsp -= lsb;
  491.         msp -= msb;
  492.     }
  493.     if (lsp > lsa)
  494.         msp--;
  495.     lsa -= lsp;
  496.     msa -= msp;
  497.     qn[i] = q;
  498.     }
  499.     /* Pack everything together again. */
  500.     msq = (qn[0] << 9) + (qn[1] >> 2);
  501.     lsq = (qn[1] << 30) + (qn[2] << 19) + (qn[3] << 8) + (qn[4] >> 3);
  502.     if (msq < 0x100000) {
  503.     msq = (msq << 1) + (lsq >> 31);
  504.     lsq <<= 1;
  505.     expt--;
  506.     }
  507.     /* We should round the quotient, but we don't. */
  508.     if (expt <= 0) {        /* Underflow.  Return 0 rather than a denorm. */
  509.     lc[msw] = sign;
  510.     lc[lsw] = 0;
  511.     } else if (expt >= 0x7ff) {    /* Overflow. */
  512.     raise(SIGFPE);
  513.     lc[msw] = sign + 0x7ff00000;
  514.     lc[lsw] = 0;
  515.     } else {
  516.     lc[msw] = sign + (expt << 20) + (msq & 0xfffff);
  517.     lc[lsw] = lsq;
  518.     }
  519.     return dc;
  520. }
  521. float
  522. __divsf3(float a, float b)
  523. {
  524.     return (float)((double)a / (double)b);
  525. }
  526.  
  527. /* ---------------- Comparisons ---------------- */
  528.  
  529. static int
  530. compared2(const long *pa, const long *pb)
  531. {
  532. #define upa ((const ulong *)pa)
  533. #define upb ((const ulong *)pb)
  534.     if (pa[msw] == pb[msw]) {
  535.     int result = (upa[lsw] < upb[lsw] ? -1 :
  536.               upa[lsw] > upb[lsw] ? 1 : 0);
  537.  
  538.     return (pa[msw] < 0 ? -result : result);
  539.     }
  540.     if ((pa[msw] & pb[msw]) < 0)
  541.     return (pa[msw] < pb[msw] ? 1 : -1);
  542.     /* We have to treat -0 and +0 specially. */
  543.     else if (!((pa[msw] | pb[msw]) & 0x7fffffff) && !(pa[lsw] | pb[lsw]))
  544.     return 0;
  545.     else
  546.     return (pa[msw] > pb[msw] ? 1 : -1);
  547. #undef upa
  548. #undef upb
  549. }
  550.  
  551. /* 0 means true */
  552. int
  553. __eqdf2(double a, double b)
  554. {
  555.     return compared2(la, lb);
  556. }
  557.  
  558. /* !=0 means true */
  559. int
  560. __nedf2(double a, double b)
  561. {
  562.     return compared2(la, lb);
  563. }
  564.  
  565. /* >0 means true */
  566. int
  567. __gtdf2(double a, double b)
  568. {
  569.     return compared2(la, lb);
  570. }
  571.  
  572. /* >=0 means true */
  573. int
  574. __gedf2(double a, double b)
  575. {
  576.     return compared2(la, lb);
  577. }
  578.  
  579. /* <0 means true */
  580. int
  581. __ltdf2(double a, double b)
  582. {
  583.     return compared2(la, lb);
  584. }
  585.  
  586. /* <=0 means true */
  587. int
  588. __ledf2(double a, double b)
  589. {
  590.     return compared2(la, lb);
  591. }
  592.  
  593. static int
  594. comparef2(long va, long vb)
  595. {                /* We have to treat -0 and +0 specially. */
  596.     if (va == vb)
  597.     return 0;
  598.     if ((va & vb) < 0)
  599.     return (va < vb ? 1 : -1);
  600.     else if (!((va | vb) & 0x7fffffff))
  601.     return 0;
  602.     else
  603.     return (va > vb ? 1 : -1);
  604. }
  605.  
  606. /* See the __xxdf2 functions above for the return values. */
  607. int
  608. __eqsf2(float a, float b)
  609. {
  610.     return comparef2(ia, ib);
  611. }
  612. int
  613. __nesf2(float a, float b)
  614. {
  615.     return comparef2(ia, ib);
  616. }
  617. int
  618. __gtsf2(float a, float b)
  619. {
  620.     return comparef2(ia, ib);
  621. }
  622. int
  623. __gesf2(float a, float b)
  624. {
  625.     return comparef2(ia, ib);
  626. }
  627. int
  628. __ltsf2(float a, float b)
  629. {
  630.     return comparef2(ia, ib);
  631. }
  632. int
  633. __lesf2(float a, float b)
  634. {
  635.     return comparef2(ia, ib);
  636. }
  637.  
  638. /* ---------------- Conversion ---------------- */
  639.  
  640. long
  641. __fixdfsi(double a)
  642. {                /* This routine does check for overflow. */
  643.     long i = (la[msw] & 0xfffff) + 0x100000;
  644.     int expt = dx(la) - dx_bias;
  645.  
  646.     if (expt < 0)
  647.     return 0;
  648.     if (expt <= 20)
  649.     i >>= 20 - expt;
  650.     else if (expt >= 31 &&
  651.          (expt > 31 || i != 0x100000 || la[msw] >= 0 ||
  652.           ula[lsw] >= 1L << 21)
  653.     ) {
  654.     raise(SIGFPE);
  655.     i = (la[msw] < 0 ? 0x80000000 : 0x7fffffff);
  656.     } else
  657.     i = (i << (expt - 20)) + (ula[lsw] >> (52 - expt));
  658.     return (la[msw] < 0 ? -i : i);
  659. }
  660.  
  661. long
  662. __fixsfsi(float a)
  663. {                /* This routine does check for overflow. */
  664.     long i = (ia & 0x7fffff) + 0x800000;
  665.     int expt = fx(ia) - fx_bias;
  666.  
  667.     if (expt < 0)
  668.     return 0;
  669.     if (expt <= 23)
  670.     i >>= 23 - expt;
  671.     else if (expt >= 31 && (expt > 31 || i != 0x800000 || ia >= 0)) {
  672.     raise(SIGFPE);
  673.     i = (ia < 0 ? 0x80000000 : 0x7fffffff);
  674.     } else
  675.     i <<= expt - 23;
  676.     return (ia < 0 ? -i : i);
  677. }
  678.  
  679. double
  680. __floatsidf(long i)
  681. {
  682.     long msc;
  683.     ulong v;
  684.     long lc[2];
  685.  
  686.     if (i > 0)
  687.     msc = 0x41e00000 - 0x100000, v = i;
  688.     else if (i < 0)
  689.     msc = 0xc1e00000 - 0x100000, v = -i;
  690.     else            /* i == 0 */
  691.     return 0;
  692.     while (v < 0x01000000)
  693.     v <<= 8, msc -= 0x00800000;
  694.     if (v < 0x10000000)
  695.     v <<= 4, msc -= 0x00400000;
  696.     while (v < 0x80000000)
  697.     v <<= 1, msc -= 0x00100000;
  698.     lc[msw] = msc + (v >> 11);
  699.     lc[lsw] = v << 21;
  700.     return dc;
  701. }
  702.  
  703. float
  704. __floatsisf(long i)
  705. {
  706.     long lc;
  707.  
  708.     if (i == 0)
  709.     lc = 0;
  710.     else {
  711.     ulong v;
  712.  
  713.     if (i < 0)
  714.         lc = 0xcf000000, v = -i;
  715.     else
  716.         lc = 0x4f000000, v = i;
  717.     while (v < 0x01000000)
  718.         v <<= 8, lc -= 0x04000000;
  719.     while (v < 0x80000000)
  720.         v <<= 1, lc -= 0x00800000;
  721.     v = ((v >> 7) + 1) >> 1;
  722.     if (v > 0xffffff)
  723.         v >>= 1, lc += 0x00800000;
  724.     lc += v & 0x7fffff;
  725.     }
  726.     return fc;
  727. }
  728.  
  729. float
  730. __truncdfsf2(double a)
  731. {                /* This routine does check for overflow, but it converts */
  732.     /* underflows to zero rather than to a denormalized number. */
  733.     long lc;
  734.  
  735.     if ((la[msw] & 0x7ff00000) < 0x38100000)
  736.     lc = la[msw] & 0x80000000;
  737.     else if ((la[msw] & 0x7ff00000) >= 0x47f00000) {
  738.     raise(SIGFPE);
  739.     lc = (la[msw] & 0x80000000) + 0x7f800000;    /* infinity */
  740.     } else {
  741.     lc = (la[msw] & 0xc0000000) +
  742.         ((la[msw] & 0x07ffffff) << 3) +
  743.         (ula[lsw] >> 29);
  744.     /* Round the mantissa.  Simply adding 1 works even if the */
  745.     /* mantissa overflows, because it increments the exponent and */
  746.     /* sets the mantissa to zero! */
  747.     if (ula[lsw] & 0x10000000)
  748.         ++lc;
  749.     }
  750.     return fc;
  751. }
  752.  
  753. double
  754. __extendsfdf2(float a)
  755. {
  756.     long lc[2];
  757.  
  758.     if (!(ia & 0x7fffffff))
  759.     lc[msw] = ia, lc[lsw] = 0;
  760.     else
  761.     extend(lc, ia);
  762.     return dc;
  763. }
  764.  
  765. /* ---------------- Test program ---------------- */
  766.  
  767. #ifdef TEST
  768.  
  769. #include <stdio.h>
  770. #include <stdlib.h>
  771.  
  772. int
  773. test(double v1)
  774. {
  775.     double v3 = v1 * 3;
  776.     double vh = v1 / 2;
  777.     double vd = v3 - vh;
  778.     double vdn = v1 - v3;
  779.  
  780.     printf("%g=1 %g=3 %g=0.5 %g=2.5 %g=-2\n", v1, v3, vh, vd, vdn);
  781.     return 0;
  782. }
  783.  
  784. float
  785. randf(void)
  786. {
  787.     int v = rand();
  788.  
  789.     v = (v << 16) ^ rand();
  790.     if (!(v & 0x7f800000))
  791.     return 0;
  792.     if ((v & 0x7f800000) == 0x7f800000)
  793.     return randf();
  794.     return *(float *)&v;
  795. }
  796.  
  797. int
  798. main(int argc, char *argv[])
  799. {
  800.     int i;
  801.  
  802.     test(1.0);
  803.     for (i = 0; i < 10; ++i) {
  804.     float a = randf(), b = randf(), r;
  805.     int c;
  806.  
  807.     switch ((rand() >> 12) & 3) {
  808.         case 0:
  809.         r = a + b;
  810.         c = '+';
  811.         break;
  812.         case 1:
  813.         r = a - b;
  814.         c = '-';
  815.         break;
  816.         case 2:
  817.         r = a * b;
  818.         c = '*';
  819.         break;
  820.         case 3:
  821.         if (b == 0)
  822.             continue;
  823.         r = a / b;
  824.         c = '/';
  825.         break;
  826.     }
  827.     printf("0x%08x %c 0x%08x = 0x%08x\n",
  828.            *(int *)&a, c, *(int *)&b, *(int *)&r);
  829.     }
  830. }
  831.  
  832. /*
  833.    Results from compiling with hardware floating point:
  834.  
  835.    1=1 3=3 0.5=0.5 2.5=2.5 -2=-2
  836.    0x6f409924 - 0x01faa67a = 0x6f409924
  837.    0x00000000 + 0x75418107 = 0x75418107
  838.    0xe32fab71 - 0xc88f7816 = 0xe32fab71
  839.    0x94809067 * 0x84ddaeee = 0x00000000
  840.    0x2b0a5b7d + 0x38f70f50 = 0x38f70f50
  841.    0xa5efcef3 - 0xc5dc1a2c = 0x45dc1a2c
  842.    0xe7124521 * 0x3f4206d2 = 0xe6ddb891
  843.    0x8d175f17 * 0x2ed2c1c0 = 0x80007c9f
  844.    0x419e7a6d / 0xf2f95a35 = 0x8e22b404
  845.    0xe0b2f48f * 0xc72793fc = 0x686a49f8
  846.  
  847.  */
  848.  
  849.  
  850. #endif
  851.